module fs_FindSolution_module_q4
  implicit none
contains
!*************************************************************************************************************!
! List of functions                                                                                           !
! real function f(x) ->set this to the function to be evaluated                                               !
! function fs_bijection( x1, x2, convergence_limit) -> finds the root of f(x) using bijection method          !
! function fs_falsePosition( x1, x2, convergence_limit) -> the root is found out by the false position method !
!**************************************************************************************************************
  real function f(x)  ! the function to be analysed
    real, intent(in)  :: x
    !print *, c
    f = sin(x) - x**3 - 1
  end function f

  recursive function fs_bijection( x1, x2, convergence_limit) result(root) !  converges to the root by the bijection method
    real, intent(in)	:: x1, x2				!  x1 and x2 should bracket the root
    real, intent(in)	:: convergence_limit	!  absolute limit for resolution of root eg.-use 0.00001
    real             :: root
    real             :: r
    
    if( f(x1) == 0 ) then
      root = x1
      !print *,'f(x1) == 0'
      return
    end if

    if ( f(x2) == 0 ) then
      root = x2
      !print *,'f(x2) == 0'
      return
    end if

    if( x1 > x2) then
       root = fs_setNan()
       !print *,'x2>x1'
       return
    end if

    r = (x2 + x1)/2.0 !  r is the possible root
    print *,r
    if ( f(r) == 0 ) then
      root = r
      !print *,'f(r) == 0'
      return
    end if

    if ( f(r) < convergence_limit .AND. f(r) > -convergence_limit) then
      root = r
      !print *,'converged'
      return
    end if

    if( f(x1)*f(r) < 0 ) then    ! if the root lies between x1 and r
      root = fs_bijection(x1, r, convergence_limit)
      return;
    end if

    if( f(r)*f(x2) < 0 ) then    ! otherwise the root lies between r and x2
      root = fs_bijection(r, x2, convergence_limit)
      return;
    end if        

    root = fs_setNan()

  end function fs_bijection

  recursive function fs_falsePosition(x1, x2, convergence_limit) result(root)
    real, intent(in)  :: x1, x2
    real, intent(in)  :: convergence_limit
    real              :: root
    real              :: y1, y2, y0
    real              :: x0
    real              :: m, c

    ! check if the points bracket the root
    if( f(x1) == 0 ) then
      root = x1
      !print *,'f(x1) == 0'
      return
    end if
    if ( f(x2) == 0 ) then
      root = x2
      !print *,'f(x2) == 0'
      return
    end if
    y1 = f(x1); y2 = f(x2);
    if( y1*y2 > 0) then
      root = fs_setNan()
      return
    end if
    if(x1>x2) then
      root = fs_setNan();
      return
    end if
    ! find the point where the line joining the two points intersect
    m = (y2 - y1)/(x2 - x1)
    c = y1 - m*x1
    x0 = -1.0*c/m
    y0 = f(x0)
    print *, x0
    ! if the root converges
    if ( y0 < convergence_limit .AND. y0 > -convergence_limit) then
      root = x0
      !print *,'converged'
      return
    end if
    ! if it does not
    if( y1*y0 < 0 ) then    ! if the root lies between x1 and r
      root = fs_falsePosition(x1, x0, convergence_limit)
      return;
    end if
    if( y0*y2 < 0 ) then    ! otherwise the root lies between r and x2
      root = fs_falsePosition(x0, x2, convergence_limit)
      return;
    end if        
    root = fs_setNan()   
  end function fs_falsePosition

  real function fs_setNan()
    real :: yo = -1
    fs_setNan = sqrt(yo)
  end function fs_setNan

end module fs_FindSolution_module_q4
